/*
 * SpecialFunctions.cs
 * 
 * Copyright (c) 2003-2005, dnAnalytics. All rights reserved.
*/

using System;
using dnAnalytics.Resources;
namespace dnAnalytics.Math
{
	///<summary>Contains special functions not found in System.Math.</summary>
	public sealed class SpecialFunctions
	{
		private SpecialFunctions() { }
		private const int MAX_ITERATIONS = 500;
		private static readonly double[] C = {
            0.99999999999999709182,
            57.156235665862923517,
            -59.597960355475491248,
            14.136097974741747174,
            -0.49191381609762019978,
            .33994649984811888699e-4,
            .46523628927048575665e-4,
            -.98374475304879564677e-4,
            .15808870322491248884e-3,
            -.21026444172410488319e-3,
            .21743961811521264320e-3,
            -.16431810653676389022e-3,
            .84418223983852743293e-4,
            -.26190838401581408670e-4,
            .36899182659531622704e-5,
        };

		private static readonly long[] FactorialTable = {
            1, 1, 2, 4, 6, 24, 120, 720, 40320, 362880, 3628800, 39916800,
            479001600, 6227020800, 87178291200, 1307674368000, 20922789888000,
            355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000
        };

		//Based on http://jakarta.apache.org/commons/math/xref/org/apache/commons/math/special/Gamma.html
		///<summary>Calculates the natural log of the gamma function evaluated at X.</summary>
		///<param name="x">value to evaluate.</param>
		///<returns>the natural log of the gamma function evaluated at X.</returns>
		///<remarks>The method returns <b>Double.NaN</b> if X &lt;= 0.</remarks>
        [System.Diagnostics.CodeAnalysis.SuppressMessage("Microsoft.Naming", "CA1706:ShortAcronymsShouldBeUppercase", MessageId = "Member")]
        public static double LnGamma(double x)
		{
			double ret;
			if (Double.IsNaN(x) || (x <= 0.0))
			{
				return Double.NaN;
			} else
			{
				double g = 607.0 / 128.0;

				double sum = 0.0;
				for (int i = 1; i < C.Length; ++i)
				{
					sum += (C[i] / (x + i));
				}
				sum += C[0];

				double tmp = x + g + .5;
				ret = ((x + .5) * System.Math.Log(tmp)) - tmp +
					  (.5 * System.Math.Log(2.0 * System.Math.PI)) + System.Math.Log(sum) - System.Math.Log(x);
			}

			return ret;

		}

		///<summary>Calculates the gamma function evaluated at X.</summary>
		///<param name="x">value to evaluate.</param>
		///<returns>the function evaluated at X.</returns>
		///<remarks>The method returns <b>Double.NaN</b> if X &lt;= 0.</remarks>
		public static double Gamma(double x)
		{
			if (Double.IsNaN(x) || (x <= 0.0))
			{
				return Double.NaN;
			} else
			{
				return System.Math.Exp(LnGamma(x));
			}
		}

		//based on Numerical Recipes 6.2
		///<summary>Calculates the incomplete gamma function evaluated at A and X.</summary>
		///<param name="a">the parameter of the gamma distribution.</param>
		///<param name="x">the integration end point.</param>
		///<returns>the incomplete gamma function evaluated at A and X.</returns>
		///<remarks>The method returns <b>Double.NaN</b> if X &lt; 0 or A &lt;= 0.</remarks>
		///<exception cref="ArithmeticException">if the the method is unable to converge to a value.</exception>
		public static double IncompleteGamma(double a, double x){
            if (Double.IsNaN(a) || Double.IsNaN(x) || (a <= 0.0) || (x < 0.0)) {
                return Double.NaN;
            } else if ( x==0 ) {
                return 0;
            } else if ( Double.IsPositiveInfinity(x) ) {
                return 1;
            } else {
                double an = 1.0/a;
                double sum = an;
                for (int n = 1; n <= MAX_ITERATIONS; ++n) {
                    an *= x / (a+n);
                    sum += an;
                    if (System.Math.Abs(an) < Double.Epsilon) {
                        return sum * System.Math.Exp(-x + a * System.Math.Log(x) - LnGamma(a));
                    }
                }
                throw new System.ArithmeticException(Strings.CalculateIncompleteGamma);
            }
        }

		///<summary>Calculates the beta function evaluated at P and Q.</summary>
		///<param name="p">value to evaluate.</param>
		///<param name="q">value to evaluate.</param>
		///<returns>the beta function evaluated at P and Q.</returns>
		///<remarks>The method returns <b>Double.NaN</b> if P,Q, P+Q &lt;= 0.</remarks>
		public static double Beta(double p, double q)
		{
			return System.Math.Exp(LnGamma(p) + LnGamma(q) - LnGamma(p + q));
		}

		///<summary>Calculates the error function evaluated at X.</summary>
		///<param name="x">value to evaluate.</param>
		///<returns>the error function evaluated at X.</returns>
		public static double Error(double x)
		{
			if (Double.IsNaN(x))
			{
				return Double.NaN;
			} else if (x == 0)
			{
				return 0;
			} else if (Double.IsPositiveInfinity(x))
			{
				return 1;
			} else if (Double.IsNegativeInfinity(x))
			{
				return -1;
			} else
			{
				return (x < 0) ? -IncompleteGamma(.5, x * x) : IncompleteGamma(.5, x * x);
			}
		}

		///<summary>Calculates the inverse error function evaluated at X.</summary>
		///<param name="x">value to evaluate.</param>
		///<returns>the error function evaluated at X.</returns>
		public static double InverseError(double x)
		{
			if (Double.IsNaN(x))
			{
				return Double.NaN;
			} else if (x == 0)
			{
				return 0;
			} else if (x == 1)
			{
				return Double.PositiveInfinity;
			} else
			{
				return (x < 0) ? -IncompleteGamma(.5, x * x) : 1 - IncompleteGamma(.5, x * x);
			}
		}

		///<summary>Calculates the factorial of the integer N.</summary>
		///<param name="n">value to evaluate.</param>
		///<returns>the factorial of the integer N.</returns>
		public static double Factorial(int n)
		{
			if (n < 0)
			{
				return Double.NaN;
			} else if (n <= 20)
			{
				return FactorialTable[n];
			} else
			{
				return System.Math.Exp(LnGamma(n + 1));
			}
		}

		///<summary>Calculates the natural log of factorial of the integer N.</summary>
		///<param name="n">value to evaluate.</param>
		///<returns>the factorial of the integer N.</returns>
		///<remarks>Returns 0 for N &lt;= 1. </remarks>
        [System.Diagnostics.CodeAnalysis.SuppressMessage("Microsoft.Naming", "CA1706:ShortAcronymsShouldBeUppercase", MessageId = "Member")]
        public static double LnFactorial(int n)
		{
			if (n < 0)
			{
				return Double.NaN;
			} else if (n <= 1)
			{
				return 0.0;
			} else if (n <= 20)
			{
				return System.Math.Log(FactorialTable[n]);
			} else
			{
				return LnGamma(n + 1);
			}
		}

		///<summary>Calculates the binomial function evaluated at N and K.</summary>
		///<param name="n">value to evaluate.</param>
		///<param name="k">value to evaluate.</param>
		///<returns>the binomial function evaluated at N and K.</returns>
		///<remarks>The method returns <b>Double.NaN</b> if N,K &lt; 0.</remarks>
		public static double Binomial(int n, int k)
		{
            checked
            {
                if (k == 0)
                {
                    return 1;
                }
                else if (k == 1)
                {
                    return n;
                }
                else if (k == n)
                {
                    return 1;
                }
                else if (k == n + 1)
                {
                    return 0;
                }
                else if (k == 2)
                {
                    return (n - 1) * n / 2.0;
                }
                else if (k == 3)
                {
                    return (n - 2) * (n - 1) * n / 6.0;
                }
                else if (k == 4)
                {
                    return (n - 3) * (n - 2) * (n - 1) * n / 24.0;
                }
                else
                {
                    return System.Math.Floor(0.5 + System.Math.Exp(LnFactorial(n) - LnFactorial(k) - LnFactorial(n - k)));
                }
            }
		}
	}
}


